|Instalar librerias necesarias (en modo silencioso -q):
[73]:
!wget -q https://raw.githubusercontent.com/arqlm/arqlm.github.io/main/_static/libraries.txt
!pip install -q -r /content/libraries.txt
!pip install -q --upgrade plotly[kaleido]
!rm -r /content/libraries.txt
!rm -r /content/sample_data/ ## Esta linea no es necesaria cuando no se trabaja en colab.google
'wget' is not recognized as an internal or external command,
operable program or batch file.
ERROR: Could not open requirements file: [Errno 2] No such file or directory: '/content/libraries.txt'
'rm' is not recognized as an internal or external command,
operable program or batch file.
'rm' is not recognized as an internal or external command,
operable program or batch file.
[74]:
## Importar librerías
import pandas as pd
import numpy as np
import base64
import datetime
import io
import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py
# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
7. Verificación de tendencias (derivas)
A veces existen tendencias sistemáticas en el valor de la variable. Por ejemplo, la concentración de un elemento puede cambiar sistemáticamente con la profundidad. En estos casos, podemos dejar de lado la idea de un dominio estrictamente estacionario utilizando la noción de cuasi-estacionariedad de segundo orden; es decir, podemos asumir que para cada ubicación, las propiedades estadísticas (y geológicas) permanecen homogéneas dentro de un vecindario.
Este vecindario puede definirse para que la suposición funcione; es decir, si existe una tendencia fuerte en los datos, se puede usar un vecindario lo suficientemente pequeño para asumir estacionariedad dentro de él.
En presencia de tendencias sistemáticas en el atributo:
Graficar medias y varianzas locales en diferentes direcciones.
Determinar si el cambio en la media o varianza es significativo en relación con la escala del modelo:
Si la tendencia no es significativa: usar los vecindarios locales para estimación y simulación, y asegurarse de revisar el modelo al final para ver qué tan bien sigue las tendencias.
Si la tendencia sigue siendo significativa dentro del vecindario, modelar la tendencia con una interpolación suave y calcular los residuales en las ubicaciones de las muestras.
Determinar si dentro de un vecindario local, de tamaño suficiente para encontrar muestras para interpolar, la tendencia sigue siendo significativa:
Actualizar las estadísticas y revisar.
Es claro que existe una relación entre la escala de la posible tendencia y el tamaño del vecindario. Además, la disponibilidad de datos de muestra es importante, ya que se necesita encontrar suficientes muestras dentro de cada pequeño vecindario para permitir la inferencia.
Como regla general, si hay abundantes muestras, el efecto de la tendencia será despreciable, ya que existirá suficiente información para informar las estadísticas locales dentro de un vecindario pequeño. Las tendencias se vuelven (muy) problemáticas cuando existen pocos datos, por ejemplo, en etapas tempranas de exploración, donde los sondajes están muy separados.
Definir dominios estacionarios es una decisión, y una muy importante, ya que determina qué datos se utilizan para inferir sobre ubicaciones no muestreadas dentro de esa región. Además, definir los vecindarios locales cuasi-estacionarios es tan importante como la definición de los dominios.
En la práctica, los dominios son dictados por el entendimiento geológico. Se realizan análisis estadísticos y espaciales para evaluar cuán apropiado es un modelo estacionario para los datos. Para este propósito se utilizan herramientas de análisis exploratorio de datos. Si los dominios son homogéneos desde el punto de vista geológico, pero muestran tendencias, esto puede manejarse definiendo el vecindario local apropiado al hacer inferencias sobre cada ubicación.
En algunos casos, esto no es suficiente y puede ser necesario modelar y remover la tendencia, con la esperanza de que los residuales obtenidos a través de este proceso sean estacionarios y capturen parte de la estructura espacial. Sin embargo, este es un equilibrio delicado y no se dará una regla general en esta etapa.
8. Cargar las bibliotecas requeridas
El siguiente código carga las bibliotecas requeridas.
[75]:
#### Instalar libreria de geoestadistica Geostatspy
!pip install geostatspy
Requirement already satisfied: geostatspy in c:\users\alvaro\anaconda3\lib\site-packages (0.0.79)
[76]:
import numpy as np # ndarrays for gridded data
import pandas as pd # DataFrames for tabular data
import geostatspy.GSLIB as GSLIB # GSLIB utilities, visualization and wrapper
import geostatspy.geostats as geostats # GSLIB methods convert to Python
import geostatspy
print('GeostatsPy version: ' + str(geostatspy.__version__))
import os # set working directory, run executables
from tqdm import tqdm # suppress the status bar
from functools import partialmethod
tqdm.__init__ = partialmethod(tqdm.__init__, disable=True)
ignore_warnings = True # ignore warnings?
import matplotlib.pyplot as plt # for plotting
import matplotlib as mpl # custom colorbar
import matplotlib.ticker as mticker # custom colorpar ticks
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
plt.rc('axes', axisbelow=True) # plot all grids below the plot elements
if ignore_warnings == True:
import warnings
warnings.filterwarnings('ignore')
import os
import matplotlib.image as mpimg
import subprocess
import time
from matplotlib.offsetbox import AnchoredText
import matplotlib.gridspec as gridspec
import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py
# trimmed inverted rainbow colormap hsv removing magenta tones
cmap = mpl.colors.ListedColormap(mpl.cm.hsv(np.linspace(0.0, 0.70, 256))).reversed()
GeostatsPy version: 0.0.79
Nuestro punto de partida sera la base de datos simplificada de la sección anterior.
[77]:
## Cargar data simplificada de la sección anterior
import pandas as pd
import numpy as np
# Cargar base de datos en pandas y aproximar a dos decimales
df = pd.read_csv('Data_sin_compositar_con_UGs.csv', sep=',', encoding='latin1').round(2).dropna()
df['UG'] = df['UG'].astype(int) # Asegurarse de que la columna 'UG' sea de tipo entero
df.head() # mostrar las primeras filas de la data
[77]:
| Este | Norte | Elevación | au_ppm | ag_ppm | cu_pct | aucn_ppm | cucn_ppm | Zmin | Alte | Lito | UG | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 472187.20 | 6925805.49 | 4213.86 | 0.30 | 2.3 | 0.01 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 |
| 2 | 472187.34 | 6925805.77 | 4211.99 | 0.47 | 16.2 | 0.03 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 |
| 3 | 472187.49 | 6925806.06 | 4210.01 | 0.31 | 2.3 | 0.02 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 |
| 4 | 472187.64 | 6925806.35 | 4208.04 | 0.29 | 2.1 | 0.01 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 |
| 5 | 472187.78 | 6925806.64 | 4206.07 | 0.40 | 2.1 | 0.01 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 |
9. Codificar base de datos en funcion de indicadores
Convendrá agregar a nuestra base de datos 3 columnas nuevas que codifiquen la variable indicador para cada UG
[78]:
# Agregar columnas de leyes para cada UG
df['ley Cu UG_1'] = np.where(df['UG'] == 1, df['cu_pct'], -99)
df['ley Cu UG_2'] = np.where(df['UG'] == 2, df['cu_pct'], -99)
df['ley Cu UG_3'] = np.where(df['UG'] == 3, df['cu_pct'], -99)
# Guardar data con leyes en formato GSLIB
GSLIB.Dataframe2GSLIB('./01_data/data_contacto.dat',df) # guardar data con leyes por UG en formato GSLIB
df.head() # mostrar las primeras filas de la data
[78]:
| Este | Norte | Elevación | au_ppm | ag_ppm | cu_pct | aucn_ppm | cucn_ppm | Zmin | Alte | Lito | UG | ley Cu UG_1 | ley Cu UG_2 | ley Cu UG_3 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 472187.20 | 6925805.49 | 4213.86 | 0.30 | 2.3 | 0.01 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 | 0.01 | -99.0 | -99.0 |
| 2 | 472187.34 | 6925805.77 | 4211.99 | 0.47 | 16.2 | 0.03 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 | 0.03 | -99.0 | -99.0 |
| 3 | 472187.49 | 6925806.06 | 4210.01 | 0.31 | 2.3 | 0.02 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 | 0.02 | -99.0 | -99.0 |
| 4 | 472187.64 | 6925806.35 | 4208.04 | 0.29 | 2.1 | 0.01 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 | 0.01 | -99.0 | -99.0 |
| 5 | 472187.78 | 6925806.64 | 4206.07 | 0.40 | 2.1 | 0.01 | -99.0 | -99.0 | OXI | ILL_CLO | IBX_MM | 1 | 0.01 | -99.0 | -99.0 |
9.1. Análisis de contacto
Definimos la siguiente función, que nos permitirá optimizar la ejecución del calculo de variogramas.
[79]:
def contacto(df, xcol, ycol, zcol, vcol1,vcol2, nlag, lagdist, lag_tol):
"""Analisis de contacto en 3D usando GSLIB GAMV."""
lag = []
npair = []
head = []
tail = []
df_ext = pd.DataFrame({"X": df[xcol], "Y": df[ycol], "Z": df[zcol],"Variable1": df[vcol1], "Variable2": df[vcol2]})
GSLIB.Dataframe2GSLIB("01_data/contacto_out.dat", df_ext)
with open("02_parametros/contacto.par", "w") as f:
f.write(" Parameters for GAMV \n")
f.write(" ******************* \n")
f.write(" \n")
f.write("START OF PARAMETERS: \n")
f.write("01_data/contacto_out.dat -file with data \n")
f.write("1 2 3 - columns for X, Y, Z coordinates \n")
f.write("2 4 5 - number of variables,col numbers \n")
f.write("0 999999 - trimming limits \n")
f.write("03_resultados/contacto.out -file for variogram output \n")
f.write(str(nlag) + " -number of lags \n")
f.write(str(lagdist) + " -lag separation distance \n")
f.write(str(lag_tol) + " -lag tolerance \n")
f.write("1 -number of directions \n")
f.write("0 90 99999 0 90 99999 -azm,atol,bandh,dip,dtol,bandv \n")
f.write("0 -standardize sills? (0=no, 1=yes) \n")
f.write("1 -number of variograms \n")
f.write("1 2 1 -tail var., head var., variogram type \n")
# Ejecutar gamv de GSLIB
result = subprocess.run(["./Gslib90/gamv", "./02_parametros/contacto.par"], capture_output=True, text=True, check=True)
with open("03_resultados/contacto.out") as f:
next(f) # skip the first line
next(f) # skip the second line
for line in f:
_, l, _, n, h,t = line.split()
lag.append(float(l))
head.append(float(h))
tail.append(float(t))
## append only if n is defined
if n.isdigit():
npair.append(float(n))
else:
npair.append(np.nan)
return lag, npair, head, tail
[80]:
lag_UG1_2, UG1_2_npair, UG1_2_head, UG1_2_tail = contacto(df,"Este","Norte","Elevación","ley Cu UG_1","ley Cu UG_2", 10, 10,10)
[81]:
from plotly.subplots import make_subplots
## Plot truncated z values against filtered average cut > 0 percentage with a line plot and marginal counter histogram
import plotly.express as px
import plotly.graph_objects as go
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
# create subplot: top = line, bottom = marginal histogram (shared x)
# enable a secondary y-axis for the bottom subplot so we can plot the bar on the right
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
vertical_spacing=0.05,
row_heights=[0.3, 0.3],
specs=[[{"type": "xy"}],
[{"type": "xy", "secondary_y": True}]])
# line trace (top)
fig.add_trace(
go.Scatter(
x=-np.array(lag_UG1_2)/2,
y=UG1_2_head,
mode='lines+markers',
name='Ley Promedio Cu [%] en UG 1',
line=dict(color='blue'),
marker=dict(size=6, color='blue')
),
row=1, col=1
)
## add z vs filtered line
# line trace (top)
fig.add_trace(
go.Scatter(
x=np.array(lag_UG1_2)/2,
y=UG1_2_tail,
mode='lines+markers',
name='Ley Promedio Cu (%) en UG 2',
line=dict(color='red'),
marker=dict(size=6, color='red')
),
row=1, col=1
)
bar = go.Bar(
x=-np.array(lag_UG1_2)/2,
y=UG1_2_npair,
name='N° celdas válidas (modelo)',
marker=dict(color='blue'),
opacity=0.3
)
bar2 = go.Bar(
x=np.array(lag_UG1_2)/2,
y=UG1_2_npair,
name='N° muestras',
marker=dict(color='red'),
opacity=0.3
)
## add bar y-axis on the right, row=2, col=1.
# add the bar to the bottom subplot using the secondary y-axis
fig.add_trace(bar, row=2, col=1)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
title_text='Celdas válidas (modelo)',
row=2, col=1,gridcolor='lightgrey'
)
fig.add_trace(bar2, row=2, col=1)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
title_text='Celdas válidas (modelo)',
row=2, col=1
)
# layout and axis labels
fig.update_xaxes(title_text='Distancia al contacto (m)', row=2, col=1)
fig.update_yaxes(title_text='Ley Promedio Cu (%)', row=1, col=1,gridcolor='lightgrey')
fig.update_yaxes(title_text='Conteo muestras', row=2, col=1)
fig.update_layout(
title='Análisis de contacto entre UG 1 y UG 2, ley de Cu (%)',
font_family="Times New Roman",
font_size=12,
font_color="black",
plot_bgcolor='white',
hovermode=False,
showlegend=True,
bargap=0.1,
height=600,
width=900
)
fig.show()
[82]:
lag_UG2_3, UGG2_3_npair, UGG2_3_head, UGG2_3_tail = contacto(df,"Este","Norte","Elevación","ley Cu UG_2","ley Cu UG_3", 10, 10,10)
[83]:
from plotly.subplots import make_subplots
## Plot truncated z values against filtered average cut > 0 percentage with a line plot and marginal counter histogram
import plotly.express as px
import plotly.graph_objects as go
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
# create subplot: top = line, bottom = marginal histogram (shared x)
# enable a secondary y-axis for the bottom subplot so we can plot the bar on the right
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
vertical_spacing=0.05,
row_heights=[0.3, 0.3],
specs=[[{"type": "xy"}],
[{"type": "xy", "secondary_y": True}]])
# line trace (top)
fig.add_trace(
go.Scatter(
x=-np.array(lag_UG2_3)/2,
y=UGG2_3_head,
mode='lines+markers',
name='Ley Promedio Cu [%] en UG 2',
line=dict(color='blue'),
marker=dict(size=6, color='blue')
),
row=1, col=1
)
## add z vs filtered line
# line trace (top)
fig.add_trace(
go.Scatter(
x=np.array(lag_UG2_3)/2,
y=UGG2_3_tail,
mode='lines+markers',
name='Ley Promedio Cu (%) en UG 3',
line=dict(color='red'),
marker=dict(size=6, color='red')
),
row=1, col=1
)
bar = go.Bar(
x=-np.array(lag_UG2_3)/2,
y=UGG2_3_npair,
name='N° muestras',
marker=dict(color='blue'),
opacity=0.3
)
bar2 = go.Bar(
x=np.array(lag_UG2_3)/2,
y=UGG2_3_npair,
name='N° muestras',
marker=dict(color='red'),
opacity=0.3
)
## add bar y-axis on the right, row=2, col=1.
# add the bar to the bottom subplot using the secondary y-axis
fig.add_trace(bar, row=2, col=1)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
title_text='Celdas válidas (modelo)',
row=2, col=1,gridcolor='lightgrey'
)
fig.add_trace(bar2, row=2, col=1)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
title_text='Celdas válidas (modelo)',
row=2, col=1
)
# layout and axis labels
fig.update_xaxes(title_text='Distancia al contacto (m)', row=2, col=1)
fig.update_yaxes(title_text='Ley Promedio Cu (%)', row=1, col=1,gridcolor='lightgrey')
fig.update_yaxes(title_text='Conteo muestras', row=2, col=1)
fig.update_layout(
title='Análisis de contacto entre UG 2 y UG 3, ley de Cu (%)',
font_family="Times New Roman",
font_size=12,
font_color="black",
plot_bgcolor='white',
hovermode=False,
showlegend=True,
bargap=0.1,
height=600,
width=900
)
fig.show()